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ABSTRACT 

We explore the possibility for the formation of Population III binaries. The 
collapse of a rotating cylinder is simulated with a three-dimensional, high- 
resolution nested grid, assuming the thermal history of primordial gas. The 
simulations are done with dimensionless units, and the results are applicable to 
low-mass as well as massive systems by scaling with the initial density. We find 
that if the initial angular momentum is as small as (3 « 0.1, where (3 is the ratio 
of centrifugal force to pressure force, then the runaway collapse of the cloud stops 
to form a rotationally-supported disk. After the accretion of the envelope, the 
disk undergoes a ring instability, eventually fragmenting into a binary. If the 
initial angular momentum is relatively large, a bar-type instability arises, result- 
ing in the collapse into a single star through rapid angular momentum transfer. 
The present results show that a significant fraction of Pop III stars are expected 
to form in binary systems, even if they are quite massive or less massive. The 
cosmological implications of Population III binaries are briefly discussed. 

Subject headings: binaries: general — cosmology: theory — hydrodynamics - 
instabilities — stars: formation 



1. Introduction 



In present-day galaxies, it is known that more than half of main-sequence stars are found 
in binaries or low-order multiple systems not only for low-mass stars but also for massive stars 
(e.g., Tohline 2002; Kroupa 2004). Regarding the upper mass limit, a very massive binary 
of the order of 1OOM is found in the Large Magellanic Cloud (Ostrov 2002). Furthermore, 
it has been revealed that the binary fraction of pre-main-sequence stars is comparable to 
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that of main-sequence stars (Mathieu 1994; Prosser et al. 1994). This implies that binaries 
possibly form prior to the pre-main-sequence phase. 

Recent studies on the formation of Population III stars have shown that the first stars 
in the universe are likely to be as massive as 100 — 1000M o (Bromm et al. 1999; Abel et al. 
2000, 2002) or to have formed in a bimodal initial mass function with peaks of ~ 1M and 
several 100M o (Nakamura & Umemura 1999, 2001). In light of present-day star formation, 
a portion of Population III stars appear able to form in binary systems. Population III 
binaries could play an important role for the thermal history of the universe through the 
formation of massive binary black holes (BHs), accreting binaries, and Population III Type 
la supernovae (SNe). For instance, the accreting binaries may emit ultraviolet radiation 
in a longer timescale than the lifetime of massive stars themselves and therefore can be a 
reionization source in the early universe (Umemura 2004). 

The binary formation for present-day stars has been studied through elaborate numerical 
simulations (Bonnell 1994; Boss 1986; Tsuribe & Inutsuka 1999; Matsumoto & Hanawa 
2003b; Nakamura & Li 2003; Machida et al. 2004). However, no study has been made on 
the formation of Population III binaries. So far, studies of Population III protostars have 
focused on the core collapse of single stars which proceeds through the cooling by primordial 
hydrogen molecules (Palla et al. 1983; Omukai & Nishi 1998; Omukai 2000). The previous 
studies imply that (1) the initial stellar core is on the order of M core = 1O~ 2 M , which is 
not so different from the present-day protostars, (2) the mass accretion rate is as high as a 
few tens of c^/G with the sound speed c s , which is numerically M ps 10 -3 — lO~ 1 M yr~ 1 
for ps 100K gas because no coolant works to cool down to less than 100K in primordial gas, 
and (3) the thermal history follows the evolution with an adiabatic exponent of 7 ~ 1.1 
until reaching the central density of n ~ 10 20 cm -3 , above which hydrogen molecules are 
dissociated and therefore the temperature increases with 7 = 5/3. 

On the basis of these findings, we study the formation of Population III binary stars. We 
perform a three-dimensional hydrodynamic simulation with nested grids in order to resolve 
the core collapse and the subsequent fragmentation. In the next section, the numerical model 
is described. In §3, we present the numerical results for the collapse and fragmentation of 
Population III stellar core. §4 is devoted to the discussion of the cosmological implications 
of Population III binaries. 
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2. Model 

In the previous three-dimensional numerical simulations of the formation of Population 
III stars (e.g., Bromm et al. 1999), filamentary structure commonly appears before it is 
broken into many pieces. Hence, we consider a cylindrical cloud as the initial condition, in 
which the cloud is assumed to be in hydrostatic balance with the central temperature and 
density of T c0 = 800K and n c0 = 1 x 10 12 cm~ 3 , respectively. The equation of state is assumed 
to be polytropic with 7 = 1.1. The polytropic constant K is proportional to n 1 ^ 1 T c0 and 
fixed to 6.87 x 10 7 in SI units throughout the calculation. Therefore, if a different initial 
density is adopted, the temperature is scaled as T c0 = 800K(n c0 /10 12 cm -3 ) 01 . We add a 
5% sinusoidal density fluctuation of one critical wavelength along the cylinder axis, above 
which the cylinder is gravitationally unstable. In addition, in the azimuthal direction, a 
perturbation of m = 2 and 4 mode is blended with Sp/p = 10%. The total mass for one 
critical wavelength is M = 1.5M for the assumed initial density. We note that, since actual 
simulations are done with dimensionless units, the present model follows the scaling as M = 
1.5M Q (n c0 /10 12 cm- 3 )( 3 T- 4 )/ 2 for the system mass and L = 1477AU(n c0 /10 12 cm- 3 )(^ 2 V 2 for 
the linear size of computational domain. For instance, the numerical results can be translated 
into a massive case with M = 423M and L = 1.04pc, if we renormalize the simulation by 
n c0 = 1 x 10 5 cm~ 3 . 

The rotation of the initial cloud is specified by a model parameter (3, which is defined by 
the ratio of the centrifugal force to the pressure gradient force, i.e., v 2 /r = — /3(l/p)(dP/dr). 
For a given /3, the specific angular momentum at the outer boundary, R, is given by j = 
[2GIR 2 P /(l + /3)] 1 / 2 , where / is the mass per unit length. We construct four models with 
different rotation of (5 — 0.0, 0.1, 0.3 and 1.0. 

The hydrodynamic evolution of the cloud is solved by a three-dimensional nested grid 
code, which is developed by Matsumoto & Hanawa (2003b). At the initial epoch, the cylinder 
is represented in (n x ,n y ,n z ) = (128, 128,32) grid points, assuming a mirror symmetry with 
respect to the equatorial plane. We introduce a new finer grid so that the resolution satisfies 
the condition that the half-width of the half-maximum density should be resolved with more 
than 8 grid points. Here 16 levels for nested grids are allowed, and therefore the final spatial 
resolution Ax = L/(128 x 2 15 ) ~ L/(A.2 x 10 6 ). Here we set L = 1477AU and therefore 
Ax = 3.5 x 10~ 4 AU. The self-gravity is calculated self-consistently by solving the Poisson 
equation with a multigrid iteration method (Matsumoto & Hanawa 2003a). 
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3. Results 

The numerical results can be divided into two cases of low and high (3 cases. Here, we 
show the typical results of these models. 

3.1. Low p case 

Figure 1 shows the density distribution and velocity fields for a case of /3 — 0.1. At 
the initial stage (Fig. la), the ratio of rotation energy E rot to gravitational energy \W\ is as 
small as -E rot /|VF| = 0.0187. The filamentary cloud fragments as a results of gravitational 
instability, and each fragment undergoes the runaway collapse. In early stages of runaway 
collapse, the dense cloud shrinks, keeping an approximately spherical shape. At the stage 
of n c = 1.04 x 10 14 cm -3 , _E rot /|IU| of spherical core reaches 0.245 because of the spin-up 
of the central cloud during the collapse. After that, the spherical collapsing core forms a 
runaway collapsing disk in the central region (Fig. 2b). When the cloud collapses near to 
the centrifugal barrier, a rotationally supported disk, the radius of which is ~ 0.035 AU is 
formed (Fig. lc) at the stage of n c = 1.82 x 10 19 cm -3 . The formation of a rotationally- 
supported disk is an important consequence for the binary formation. If the equation of 
state is perfectly isothermal (7 = 1), the runaway collapse of a rotating cloud continues to 
form the central singularity without forming a rotationally supported disk (Saigo & Hanawa 
1998; Saigo et al. 2000). However, the present simulation shows that the runaway collapse 
stops before producing the central singularity in the case of 7 = 1.1. This is an essential 
difference between isothermal and 7 = 1.1 evolution. After a rotationally supported disk 
forms, the envelope accretes onto the disk (Fig. Id). In this disk, the Toomre Q(= K,c s /nGa) 
is slightly less than unity. However, the disk radius is smaller than the critical wavelength 
of ring instability. The disk is therefore stable at this stage. The disk radius and mass grow 
via accretion from the infalling envelope. 

When the radius increases up to 0.1 AU, which is larger than the critical wavelength, 
the disk suffers from ring instability (Fig. le). The rotating ring is unstable against nonax- 
isymmetric modes, finally evolving into a binary core (Fig. If). The binary core separation is 
0.1 AU. The simulation was terminated just before the central density reaches ~ 10 20 cm -3 . 

The subsequent equation of state is expected to deviate from 7 = 1.1 to 5/3. When the 
binary core forms, the core mass is approximately 1O~ 2 M . The binary core can continue 
to grow through further accretion. The accretion rate onto the binary system is M ps 
1 x lO -1 M yr -1 . Saigo & Hanawa (1998) derived a Larson- Penston type accretion rate onto 
a disk as M = 8.15c 3 /G. The present accretion rate is comparable to this prediction for the 
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temperature of ~ 5,000K which is attained by the stage of n ~ 1 x 10 17 cm -3 . Except for 
the formation of a binary core, the mass accretion after core formation is likely to be similar 
to the case of a single Population III star (Omukai & Nishi 1998). 

For another low (3 case, we have simulated a (3 = 0.3 model as well. As a result, we 
have found that the dynamical evolution is very similar to the case of (3 — 0.1. Eventually, 
a binary core forms. 

3.2. High (3 case 

Models with (3 = 1.0 exhibit the evolution different from the low (3 models described 
in the previous section. Figure 2 shows the evolution of a high (3 case. At the initial 
stage (Fig. 2a), £ , ro t/|VF| = 0.103. The runaway collapse proceeds in a manner similar to 
the case of f3 — 0.1 (Fig. 2b). Until a nearly spherical collapsing cloud forms, -E rot /|W^| 
increases up to -E rot /|VF| = 0.327, owing to the enhanced rotation energy. Then a bar 
mode structure emerges as shown in Figure 2c, because the bar mode instability can set in 
when i? rot /|VF| > 0.27 (see e.g., Pickett et al. 1996; Toman et al. 1998). Simultaneously, the 
angular momentum is transferred by a bar mode, as studied by Imamura, Durisen, & Pickett 
(2000) and Saigo et al. (2002). The ambient matter accretes onto this bar (Fig. 2d), but no 
ring forms owing to effective angular momentum transfer. As a result, the bar collapses into 
a denser bar (Fig. 2e), eventually reaching stellar density of 10 20 cm~ 3 (Fig. 2f) without 
fragmenting into a multiple system. 

To see the dependence on the initial perturbation, we ran simulations with Sp/p = 
0, 10%, and 50%. As a result, we found that the final results are almost the same regardless 
of the initial perturbations. This is because the initial azimuthal perturbation is smoothed 
out in the course of the thermal history with 7 = 1.1. 

4. Conclusions and Discussion 

We have simulated the collapse of a rotating cylinder while assuming the primordial gas 
thermal history. As a result, we have found that to form a binary system, the formation 
of a rotationally supported disk is indispensable, which is realized only by a lower rotation 
parameter. The rotation-supported disk undergoes ring-mode instability, eventually frag- 
menting in to binary cores. In contrast, the case with a relatively high rotation parameter 
leads to a bar-mode instability and the angular momentum efficiently transferred by the bar, 
resulting in the collapse into a single stellar core. These results show that the initial condi- 
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tions with less angular momentum are more prone to binary formation and therefore imply 
that a significant fraction of Population III stars are expected to form in binary systems. 

In the present simulation, we have started with high initial density. If we use the scaling 
for mass and linear size and renormalize the physical quantities by a lower initial density 
such as n c0 < 1 x 10 5 cm~ 3 , then the simulation predicts the formation of a massive binary 
with several hundred M & and a separation of ~ 100AU. Nakamura & Umemura (2001) 
showed that the mass of Population III stars could bifurcate according to the initial filament 
density. If the initial density is above the critical density of H2 molecules, n C Q > 10 5 cm~ 3 , 
then the fragment mass is peaked around 1M , whereas the peak is shifted to pa 100M o for 
n c0 < 1 x 10 5 cm" 3 . 

The present simulations show that the binary formation can occur in the wide mass 
range of 1 to several xlOOM , if we start with the initial density of 10 3 - 10 12 cm -3 . 

The recent theoretical analyses of the evolution of metal-free stars predict that the fate 
of Population III stars can be classified as follows (e.g., Heger & Woosley 2002): 

1. A star of 1M Q < M < 8M Q evolves into a white dwarf after asymptotic giant branch 
mass loss. 

2. A star of 8M < M < 25M results in a Type II SN. 

item A star of 25M < M < 14OM probably collapses into a BH. 

3. A star of 14OM < M < 26OM is partly or completely disrupted by the electron- 
positron pair instability. 

4. A star with mass of M > 26OM collapses completely to a BH without ejecting any 
heavy elements. 

Combining such Population III stellar evolution with the present results, we may expect 
a wide variety of Population III binaries. If a binary is as massive as > 3OOM , it can 
result in a binary BH. Recently, Belczynski et al. (2004) argued that a Population III binary 
BH could be the intensive gravitational wave burst sources, which may be detected by the 
Laser Interferometer Gravitational- Wave Observatory. If there is a slight mass difference in 
a massive binary, a system composed of a ~ 1OOM Population III star and several 1OOM 
BHs may form. Then the mass ejected from a massive star can accrete onto a companion 
BH, resulting in a possible source of ionizing radiation in the early universe (Umemura 2004). 
The Population III BH accretion may be an important clue for the early reionization that is 
suggested by the Wilkinson Microwave Anisotropy Probe (Ricotti & Ostriker 2004a,b). If a 
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low-mass Population III binary forms, Population III Type la SN might occur. In this case, 
heavy elements synthesized in a more rapidly evolving star can pollute a companion star by 
stellar wind. The evolution is likely to be sensitive to the metallicity-dependent opacity of 
accretion flow. Population III Type la SNe seem worth studying in the future, because they 
may play an important role in the thermal history in the universe. 
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Fig. 1. — Density and velocity distributions in the x-y plane on the midplane (bottom panels) 
and x-z plane of y = (top panels) for a low (3 model with (3 = 0.1 at six different stages. 
The color scale and contour curves denote the density in the logarithmic scale. Arrows 
show the velocity distributions. The maximum density and time, (n max [cm -3 ], t [yr]) is 
(a)(1.05xl0 12 , 0.0), O)(8.82xl0 17 , 2.1537xl0 2 ), (c)(1.82xl0 19 , 2.1545xl0 2 ), (d)(6.80x 10 19 , 
2.1549 x 10 2 ), (e)(2.35 x 10 19 , 2.1557 x 10 2 ), and (J) (1.32 x 10 21 , 2.1574 x 10 2 ) at each stage. 
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Fig. 2. — Same as Fig. 1 but for a high /3 model with (3 = 1.0. The maximum density and 
time, (n max [cm" 3 ], t [yr]) is (a)(1.05 x 10 12 , 0.0), (6)(9.93 x 10 15 , 1.8346 xlO 2 ), (c)(1.05xl0 
1.8407 xlO 2 ), (<i)(6.89xl0 17 , 1.8453 x 10 2 ), (e)(4.42 x 10 18 , 1.8483 x 10 2 ), and (/)(6.66xl0 
1.8492 x 10 2 ) at each stage. 
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